
################# ------------------------------------------------------------------------------------- #################
################# -------------------------- *** COMPUTATION OF MOMENTS *** --------------------------- #################
################# ------------------------------------------------------------------------------------- #################


#---------------------------------------------------------------------------------------------------------
#FUNCTION TO COMPUTE UNCONDITIONAL MOMENTS (WINDOWS OUTSIDE A DEFAULT) & SCATTER PLOTS BASED ON SIMULATED DATA
function compute_moments(ce::CE_Economy_M, y_sim, type_sim, b_sim, π_sim, def_sim, def_pol, def_id, q_sim, SP_sim, SP_HR_sim, c_sim, tb_sim, m_sim, ζ_sim, dlnSP_sim, dBE_sim, b_sim_tmp, inf_bp; do_plots=true, WINDOWS=1)
#Note: Inputs of this function should have the same ordering as the output from function "simul_econ"
#---------------------------------------------------------------------------------------------------------

#Reputation premium
RP_sim  = SP_sim-SP_HR_sim

#----------------------------------------------------------------------------------------------------
# Filter only by default episodes
#----------------------------------------------------------------------------------------------------

        def_fre  =  mean(def_pol)              *400  #pps
        SP_mean  =  mean(SP_sim[def_sim.==0])  *100  #bps; already annualized
        SP_std   =  std(SP_sim[def_sim.==0])   *100  #bps; already annualized
        by       =  b_sim./y_sim
        by_mean  =  mean(by[def_sim.==0])      *100  #pp
        Bv_mean  =  mean(ce.Bval ./ (ce.Bval+mean(b_sim[def_sim.==0])*((1-ce.mb)*ce.zb + ce.mb)) )

#----------------------------------------------------------------------------------------------------
# Exclude periods in which: Government is in default & +-Q quarters Window
#----------------------------------------------------------------------------------------------------
if WINDOWS==1
        #Choose window size (quarters)
        window_l  = 8    #Drop observations (t) in which Default(t +- window)==1
        window_r  = 4

        #Filter sims by no-default episodes
        tmp_def   = def_id[def_sim.==0]
        tmp_y     = log.(y_sim[def_sim.==0])
        tmp_c     = log.(c_sim[def_sim.==0])
        tmp_tb    = tb_sim[def_sim.==0]./(y_sim[def_sim.==0])
        tmp_SP    = SP_sim[def_sim.==0]
        tmp_SP_zH = SP_HR_sim[def_sim.==0]
        tmp_RP    = RP_sim[def_sim.==0]
        tmp_π     = π_sim[def_sim.==0.]
        tmp_type  = type_sim[def_sim.==0.]

        # Maximum number of windows
        Nwnd = round(Int64, maximum(tmp_def))

        ### Create matrices to fill in
        std_SP           = zeros(Nwnd).*NaN
        #Table 6
        std_c_gdp        = zeros(Nwnd).*NaN
        std_tb_gdp       = zeros(Nwnd).*NaN
        cor_c_gdp        = zeros(Nwnd).*NaN
        cor_tb_gdp       = zeros(Nwnd).*NaN
        cor_SP_gdp       = zeros(Nwnd).*NaN
        #Table 7 [Panel A]
        mean_π           = zeros(Nwnd).*NaN
        std_π            = zeros(Nwnd).*NaN
        cor_π_gdp        = zeros(Nwnd).*NaN
        #Table 8
        mean_RP          = zeros(Nwnd).*NaN
        mean_RPSP        = zeros(Nwnd).*NaN
        std_RP           = zeros(Nwnd).*NaN
        std_SP_zH        = zeros(Nwnd).*NaN
        mean_RPSP_Low_y  = zeros(Nwnd).*NaN
        cor_RP_gdp       = zeros(Nwnd).*NaN
        cor_RPSP_gdp     = zeros(Nwnd).*NaN
        #Filters
        y_low            = (mean(y_sim[def_sim.==0]) - std(y_sim[def_sim.==0]))
        y_high           = (mean(y_sim[def_sim.==0]) + std(y_sim[def_sim.==0]))
        U_out, Y_out, SP_out = [], [], [];


        #Compute statistics in no-default windows (+-2 years) [by window]
        #---------------------------------------------------------------
        for t=1:Nwnd

        #Construct windows [X,Y] days after/before a default
        ind_first  = searchsortedfirst(def_id[def_sim.==0],t-1)   # First day after a default
        ind_last   = searchsortedfirst(def_id[def_sim.==0],t)     # Last day before a default
        period_tmp = ind_first+window_l:ind_last-window_r         # Keep only days before/after a Y window.

        if length(period_tmp)>20
                std_SP[t]      = std( tmp_SP[period_tmp])
        #Table 6 - standard business cycles moments)
                cor_SP_gdp[t] = cor(tmp_SP[period_tmp], tmp_y[period_tmp]  )
                std_c_gdp[t]  = std(tmp_c[period_tmp])  ./ std(tmp_y[period_tmp]  )
                std_tb_gdp[t] = std(tmp_tb[period_tmp]) ./ std(tmp_y[period_tmp]  )
                cor_c_gdp[t]  = cor(tmp_c[period_tmp],  tmp_y[period_tmp]  )
                cor_tb_gdp[t] = cor(tmp_tb[period_tmp], tmp_y[period_tmp]  )

        #Filters; high and low output
                High_y            = exp.(tmp_y[period_tmp]).>=y_high
                High_y            = ifelse.(High_y.==0, NaN, High_y)
                Low_y             = exp.(tmp_y[period_tmp]).< y_low
                Low_y             = ifelse.(Low_y.==0, NaN, Low_y)

         #Table 7 - Misreports, BE, SP [quarterly frequency]
                type_I            = ifelse.(tmp_type[period_tmp].==1, NaN, 1)
                mean_π[t]         = mean( filter(!isnan, tmp_π[period_tmp]     .* type_I) )       #Quarterly rate
                std_π[t]          = std(  filter(!isnan, tmp_π[period_tmp]     .* type_I) )       #Quarterly rate
                if isnan.(mean_π[t])==false && ce.nz>2
                cor_π_gdp[t]      = cor(  filter(!isnan, tmp_π[period_tmp].* type_I), filter(!isnan, tmp_y[period_tmp].* type_I) )
                end
 
       #Table 8  - Reputation premium
                mean_RP[t]         = mean(filter(!isnan, tmp_RP[period_tmp] ))
                mean_RPSP[t]       = 100*mean(filter(!isnan, tmp_RP[period_tmp]./tmp_SP[period_tmp] ))            #Share of spread explained by RP
                std_RP[t]          = std(filter(!isnan,  tmp_RP[period_tmp] ))      ./ std_SP[t]
                std_SP_zH[t]       = std(filter(!isnan,  tmp_SP_zH[period_tmp] ))   ./ std_SP[t]
                mean_RPSP_Low_y[t] = 100*mean(filter(!isnan, tmp_RP[period_tmp]./tmp_SP[period_tmp] .*Low_y ))    #Share of spread explained by RP, Low Output
                cor_RP_gdp[t]      = cor(tmp_RP[period_tmp], tmp_y[period_tmp]  )
                cor_RPSP_gdp[t]    = cor(tmp_RP[period_tmp]./tmp_SP[period_tmp], tmp_y[period_tmp])

                #Simulated Data for the Scatter Plot
                if do_plots==true
                    append!(Y_out, exp.(tmp_y[period_tmp]))                   #Output
                    append!(SP_out, tmp_SP[period_tmp])                       #Spread
                    append!(U_out,  tmp_RP[period_tmp]./tmp_SP[period_tmp])   #Reputation Premium / Spread
                end

        end
        end

        #------------------------------------------------------
        #Compute means across windows
        #------------------------------------------------------
        #Table 6: Untargeted Moments - Standard business cycles moments - 
        _std_c_gdp  = mean(std_c_gdp[ .!isnan.(std_c_gdp) ])
        _std_tb_gdp = mean(std_tb_gdp[.!isnan.(std_tb_gdp)])
        _cor_c_gdp  = mean(cor_c_gdp[ .!isnan.(cor_c_gdp) ])
        _cor_tb_gdp = mean(cor_tb_gdp[.!isnan.(cor_tb_gdp)])
        _cor_SP_gdp = mean(cor_SP_gdp[.!isnan.(cor_SP_gdp)])

        #Table 7 [panel A]: Untargeted Moements - Misreports, BE, Spreads
        _mean_π        = mean(mean_π[     .!isnan.(mean_π)    ])*100   #Quarterly rate
        _std_π         = mean(std_π[      .!isnan.(std_π)     ])*100   #Quarterly rate
        _cor_π_gdp     = mean(cor_π_gdp[  .!isnan.(cor_π_gdp)])

        #Table 8 - The reputation premium
        _mean_RP         = mean(mean_RP[       .!isnan.(mean_RP)])            #E(RP)
        _mean_RPSP       = mean(mean_RPSP[     .!isnan.(mean_RPSP)])          #E(RP/SP)
        _std_RP          = mean(std_RP[        .!isnan.(std_RP)])             #σ(RP)/σ(SP)
        _std_SP_zH       = mean(std_SP_zH[     .!isnan.(std_RP)])             #σ(RP)/σ(SP) ∣ ζH
        _mean_RPSP_Low_y = mean(mean_RPSP_Low_y[ .!isnan.(mean_RPSP_Low_y)])  #E(RP/SP) ∣ Low Y
        _cor_RP_gdp      = mean(cor_RP_gdp[.!isnan.(cor_RP_gdp)])             #cor(RP,logY)
        _cor_RPSP_gdp    = mean(cor_RPSP_gdp[.!isnan.(cor_RPSP_gdp)])         #cor(RP/SP,logY)

        if do_plots==true
            U_out = convert(Array{Float64,1}, U_out)
            Y_out = convert(Array{Float64,1}, Y_out)
            SP_out= convert(Array{Float64,1}, SP_out)

            #SCATTER: RP vs log Output
            p3 = scatter(log.(Y_out), U_out, ylim=(0., 0.7), xlabel="log(y)", ylabel="Reputation premium over spreads", label="", legend=:topright, size=(600, 600), dpi=150, xtickfont=font(15), ytickfont=font(15), legendfont=font(15))
            savefig(p3, "figures/simulations/scatter_Output_RPSP.png")


            #SCATTER: RP vs Spreads
            p4 = scatter((SP_out), U_out, ylim=(0., 0.7), xlim=(3.,16.), xlabel="Spreads", ylabel="Reputation premium over spreads", label="", legend=:topright, size=(600, 600), dpi=150, xtickfont=font(15), ytickfont=font(15), legendfont=font(15))
            savefig(p4, "figures/simulations/scatter_SP_RPSP.png")

        end
end

#Returns stats
if WINDOWS==1 #default option
return  def_fre,SP_mean,  SP_std,  by_mean, Bv_mean,                #Table 5 - Targeted moments
        _std_c_gdp,_std_tb_gdp,_cor_c_gdp, _cor_tb_gdp,_cor_SP_gdp, #Table 6
        _mean_π,_std_π, _cor_π_gdp,                                  #Table 7 [panel A]
        _mean_RP,_mean_RPSP,_std_RP,_std_SP_zH,_mean_RPSP_Low_y, _cor_RP_gdp,_cor_RPSP_gdp #Table 8
       
else #case in which we don't compute moments around the windows
return def_fre,SP_mean,SP_std,by_mean,Bv_mean
end
end
#---------------------------------------------------------------------------------------------------------
